Predict with sims from density and diet models to propagate uncertainty
# Predict and then for loop to summarise and avoid having too large pred grids in the memorynsim <-500sim_dens <-predict(mcod, newdata = pred_grid, nsim = nsim) |>as.data.frame()sim_her <-predict(mher, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))sim_spr <-predict(mspr, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))sim_sad <-predict(msad, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))
The density prediction is used later for spatial overlap + it gets left_joined to the other sims, so it gets a special treatment
pred_dens_space <- pred_grid_density |>filter(year %in%c(1994, 2022)) |>summarise(density_mean =mean(sim_density),.by =c(year, X, Y) )annotations <-data.frame(xpos =c(-Inf, -Inf),ypos =c(-Inf, Inf),year =c(1994, 2022),hjust =c(-0.3, -0.3),vjust =c(-29.2, 1.5))# Max density by species (because I trim the plot)pred_dens_space |>summarise(max =max(density_mean))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
ggsave(paste0(home, "/figures/spatial_cod.pdf"), width =17, height =9, units ="cm")
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Plot annual predation with uncertainty
# Summarise predator biomass by year and sim and join that to calculate per capita predationpred_dens <- pred_grid_density |>summarise(cod_biomass =sum(sim_density *9), # area is 9 km2.by =c(year, sim) )# Get the predation by year in a loop, else it becomes too long with all the sims, running out of memory...herlist <-list()sadlist <-list()sprlist <-list()tic()for (i inunique(pred_grid$year)) { sim_her_sub <- sim_her |>filter(year == i) sim_sad_sub <- sim_sad |>filter(year == i) sim_spr_sub <- sim_spr |>filter(year == i) pred_grid_density_sub <- pred_grid_density |>filter(year == i) herlist[[i]] <-get_annual_predation(sim_her_sub, prey_name ="Herring") sadlist[[i]] <-get_annual_predation(sim_sad_sub, prey_name ="Saduria") sprlist[[i]] <-get_annual_predation(sim_spr_sub, prey_name ="Sprat")gc()}tot_pred <-bind_rows(bind_rows(herlist),bind_rows(sadlist),bind_rows(sprlist))toc()
275.448 sec elapsed
# Prepare data for fitting gams and for plottingclean_pred <- tot_pred |>mutate(species =as.factor(species))# Fit gamma-GAMs to annual predation metricsm_pop <-brm( pred_median ~ species +s(year, by = species),family =Gamma(link ="log"),control =list(adapt_delta =0.99),data = clean_pred |>mutate(pred_median = pred_median /1000))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
p1 + p2 + p3 +plot_layout(axes ="collect")
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
ggsave(paste0(home, "/figures/spatial_predation.pdf"), width =19, height =12, units ="cm")
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Plot CV in space for cod density and predation
pred_grid_density_sub <- pred_grid_density |>filter(year %in%c(1994, 2022))tic()cv_pred <-bind_rows(get_cv_predation(sim_her, years =c(1994, 2022), prey_name ="Herring"),get_cv_predation(sim_sad, years =c(1994, 2022), prey_name ="Saduria"),get_cv_predation(sim_spr, years =c(1994, 2022), prey_name ="Sprat"))toc()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
tag_facet(p1, fontface =1, col ="gray30")
Warning: Removed 1136 rows containing missing values or values outside the scale range
(`geom_raster()`).
ggsave(paste0(home, "/figures/supp/cv_density_relative_prey.pdf"), width =19, height =10, units ="cm")
Warning: Removed 1136 rows containing missing values or values outside the scale range
(`geom_raster()`).
Plot relative prey weight predictions in space for two years
# Here we can use the get_cv_predation and just use the median# Max predation by species (if I trim the plot)cv_pred |>summarise(max =max(mean), .by = species)
# A tibble: 3 × 2
species max
<chr> <dbl>
1 Herring 0.00687
2 Saduria 0.0272
3 Sprat 0.0444
# Plot!viridis(n =3, option ="magma") # high color
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
p1 + p2 + p3 +plot_layout(axes ="collect")
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
ggsave(paste0(home, "/figures/spatial_relative_prey.pdf"), width =19, height =12, units ="cm")
Warning: Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Removed 284 rows containing missing values or values outside the scale range
(`geom_raster()`).
Spatial overlap over time
# Calculate overlap indices in spacepred_grid_sum <- pred_grid_density |>mutate(cod_spr_ovr =loc_collocspfn(pred = sim_density, prey = biomass_spr),cod_her_ovr =loc_collocspfn(pred = sim_density, prey = biomass_her),cod_sad_ovr =loc_collocspfn(pred = sim_density, prey = saduria),.by =c(year, sim) ) |># Stop here for plotting overlap in spacesummarise(cod_spr_ovr_tot =sum(cod_spr_ovr),cod_her_ovr_tot =sum(cod_her_ovr),cod_sad_ovr_tot =sum(cod_sad_ovr),.by =c(year, sim) ) |># Now calculate quantiles of annual overlap indices across simssummarise(cod_spr_ovr_tot_median =quantile(cod_spr_ovr_tot, prob =0.5),cod_spr_ovr_tot_lwr =quantile(cod_spr_ovr_tot, prob =0.1),cod_spr_ovr_tot_upr =quantile(cod_spr_ovr_tot, prob =0.9),cod_her_ovr_tot_median =quantile(cod_her_ovr_tot, prob =0.5),cod_her_ovr_tot_lwr =quantile(cod_her_ovr_tot, prob =0.1),cod_her_ovr_tot_upr =quantile(cod_her_ovr_tot, prob =0.9),cod_sad_ovr_tot_median =quantile(cod_sad_ovr_tot, prob =0.5),cod_sad_ovr_tot_lwr =quantile(cod_sad_ovr_tot, prob =0.1),cod_sad_ovr_tot_upr =quantile(cod_sad_ovr_tot, prob =0.9),.by = year )clean_ovr <- pred_grid_sum |>pivot_longer(c( cod_spr_ovr_tot_median, cod_spr_ovr_tot_lwr, cod_spr_ovr_tot_upr, cod_her_ovr_tot_median, cod_her_ovr_tot_lwr, cod_her_ovr_tot_upr, cod_sad_ovr_tot_median, cod_sad_ovr_tot_lwr, cod_sad_ovr_tot_upr ),names_to ="group", values_to ="value" ) |>mutate(species =ifelse(grepl("her", group), "Herring", "Sprat"),species =ifelse(grepl("sad", group), "Saduria", species),var =ifelse(grepl("median", group), "median", "error") ) |># dplyr::select(-group) |># pivot_wider(values_from = value, names_from = var) |>mutate(species =as.factor(species))# Fit Beta gams to time series of overlapm <-brm( value ~ species +s(year, by = species),family =Beta(),control =list(adapt_delta =0.99),data = clean_ovr |>filter(var =="median"))